from mechine import Heat_recovery_system, steam_turbine, Condenser

Ne_DT = 346.625519
Pt_DT = 15.479736
Dsup_DT = 393.6354572222222
Sup_wat_T = 536.259888
Pb_DT = 5.431763
Re_attw_D = 0
Pre_real = 2.584305
Re_wat_T = 517.691833
h = 0.2
n_opt_initial = 40.050761083333335
Cond_wat_T = 51.08432
Cond_wat_D = 307.31041805555554
Feed_wat_D = 411.499905
Ta = 3.606827

import numpy as np
import os
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei']  # 指定默认字体：解决plot不能显示中文问题
plt.rcParams['axes.unicode_minus'] = False  # 解决保存图像时负号'-'显示为方块的问题

# 计算净功率
def Net_power(Ne, n_opt):
    power_consumption_fan = (n_opt / 110) ** 3 * 115 * 0.001 * 56
    current_net_power = Ne - power_consumption_fan
    return power_consumption_fan, current_net_power

def sim_run(simtime, h=0.01, sample_time=5, state_result=None, ita_current=0.9, kc_current=0.04, list=[], nelist=[]):
    global Pb_DT
    if state_result is None:
        state_result = []

    turbine = steam_turbine(ita_current)
    heatsystem = Heat_recovery_system()
    coldsystem = Condenser()

    Ne, Ps1, Ps2, P1, P2, T1, T2, Tj, Ty2, Ps3, \
        Ps4, Ps5, Ps6, Ps7, Tn7, Tl7, Tw7, hd7, Dd6, hd6, Tn6, \
        Tl6, Tw6, Dd5, hd5, Tn5, Tl5, Tw5, Dd3, hd3, Tn3, Tl3, Tw3, \
        Dd2, hd2, Tn2, Tl2, Tw2, Dd1, hd1, Tn1, Tl1, Tw1, Twc_m, Tw4, \
        drhol4_dt, rhol4_m, Pb = (state_result[item] for item in range(len(state_result)))

    # 自找数据
    # Ne_DT = data[1]  # 机组负荷
    # Pt_DT = data[2]  # 主汽压力
    # Dsup_DT = data[3]  # 主汽流量
    # Pre_real = data[4]  # 再热蒸汽压力
    # Re_wat_T = data[5]  # 再热蒸汽温度
    # Cond_wat_T = data[6]  # 凝结水温度
    # Cond_wat_D = data[7]  # 凝结水流量
    # Feed_wat_D = data[8]  # 主给水流量
    # Sup_wat_T = data[9]  # 主汽温度
    # Ta = data[10]  # 环境温度
    # # n_num = data[11]  # 风机运行台数
    # n_opt_initial = data[11]  # 风机转速, 二分法数据没有风机台数
    # # Pb_DT = data[-1] / 1000  # 背压
    # Re_attw_D = 0
    # # if Pb != 0.015:
    # #     Pb_DT = Pb
    #Pb_DT = Pb
    print(f'pbdt:{Pb_DT}')
    print(f'pb:{Pb}')

    for j in range(int(sample_time / h)):
        # 汽轮机仿真
        Ds, pt, ht, D1, p1, D1j, h1, D2, p2, Dj2, h2, P2, T2, D2r, D3, p3, h3, Dj3, D4, p4, h4, \
            Dj4, D5, p5, h5, Dj5, D6, p6, h6, Dj6, D7, p7, h7, Dj7, hc, Dc, P1, T1, Tj, Ty2, Ne, Cond_T, Ne8 = turbine.simulation(
            Sup_wat_T, Pt_DT, Dsup_DT, Ps1, Ps2, Re_attw_D,
            Ne_DT, P1, P2, Pre_real, T1, T2, Re_wat_T, Tj, Ty2,
            Ps3, Ps4, Ps5, Ps6, Ps7, Pb, h)

        # 凝汽器仿真
        _, Pb, Dec, Tc, Hc1, Vf, NTU, Qc = coldsystem.simulation(hc, Dc, Pb, n_opt_initial, Ta, h, kc_current)

        # 回热系统仿真
        Ps7, Dd7, hd7, Tn7, Tl7, Tw7, Ps6, Dd6, hd6, Tn6, Tl6, Tw6, Ps5, Dd5, hd5, Tn5, Tl5, Tw5,\
        Ps4, Tw4, Ps3, Dd3, hd3, Tn3, Tl3, Tw3, Ps2, Dd2, hd2, Tn2, Tl2, Tw2, Ps1, Dd1, hd1, Tn1, Tl1, Tw1,\
        rhol4, Twc_m, Tw4_m, rhol4_m, drhol4_dt = heatsystem.simulation(
            Cond_wat_T, Tn7, Tl7, Ps7, Tw7, D7, h7, hd7, Dd6, hd6,Tn6, Tl6, Ps6, Tw6, D6, h6, Dd5, hd5,Tn5, Tl5, Ps5,
            Tw5, D5, h5, Ps4, Dd3, hd3, D4, h4, Dsup_DT, Tn3, Tl3, Ps3, Tw3, D3, h3, Dd2, hd2,Tn2, Tl2, Ps2, Tw2,
            D2, h2, Dd1, hd1, Tn1, Tl1, Ps1, Tw1, D1, h1,Cond_wat_D, Feed_wat_D, Tw4,Twc_m,drhol4_dt,rhol4_m, h)
        list.append(Pb)
        nelist.append(Ne)

    state_result = [Ne, Ps1, Ps2, P1, P2, T1, T2, Tj, Ty2, Ps3,
                    Ps4, Ps5, Ps6, Ps7, Tn7, Tl7, Tw7, hd7, Dd6, hd6,
                    Tn6, Tl6, Tw6, Dd5, hd5, Tn5, Tl5, Tw5, Dd3, hd3,
                    Tn3, Tl3, Tw3, Dd2, hd2, Tn2, Tl2, Tw2, Dd1, hd1,
                    Tn1, Tl1, Tw1, Twc_m, Tw4, drhol4_dt, rhol4_m, Pb]

    mid_vars = [Vf, NTU, Qc, Dec]


    return state_result, mid_vars, list, nelist

def optimize_for_speed(m1, simtime, h, sample_time, tolerance, ita_current, kc_current, stable_state=None):
    global Ne_DT
    global Pb_DT
    if stable_state is None:
        # state_result = [
        #     330, 3.571, 2.261, 2.1, 2, 480, 500, 520, 900, 1.203,
        #     0.6006, 0.3668, 0.1441, 0.04862, 62, 70, 77.85, 253, 5, 350,
        #     83, 95, 107, 5, 474, 113, 120, 138, 5, 712,
        #     167, 180, 188, 5, 823, 200, 210, 218, 5, 963,
        #     220, 236, 245, 54, 180, 0.01, 0.51, 0.00975
        # ]
        state_result = [Ne_DT, 3.571, 2.261, 2.1, 2, 480, 500, 520, 900, 1.203,
                        0.6006, 0.3668, 0.1441, 0.04862, 62, 70, 77.85, 253, 5, 350,
                        83, 95, 107, 5, 474, 113, 120, 138, 5, 712,
                        167, 180, 188, 5, 823, 200, 210, 218, 5, 963,
                        220, 236, 245, 54, 180, 0.01, 0.51, Pb_DT / 1000]
        pb_list = [Pb_DT / 1000]
        print(f'初始状态:{state_result}')
    else:
        state_result = stable_state
        pb_list = [state_result[-1]]
        Ne_DT = stable_state[0]
        Pb_DT = state_result[-1]
    # print(f'输入数据：{data1}')

    ne_list = [Ne_DT]
    while True:
        state_result, mid_vars, pb_list, ne_list = sim_run(simtime, h=h, sample_time=sample_time, state_result=state_result, ita_current=ita_current, kc_current=kc_current, list=pb_list, nelist=ne_list)  # simtime 25
        result_Pc_m1 = state_result[-1] * 1000  # 计算背压
        print(f'每步状态{state_result}')

        # 判断误差是否满足条件
        if abs((result_Pc_m1 - Pb_DT) / Pb_DT) < tolerance:
            Ne_m1 = state_result[0]  # 负荷
            power_consumption_fan_m1, net_power_m1 = Net_power(float(Ne_m1), m1)
            break  # 满足误差条件，退出循环
        else:
            Pb_DT = result_Pc_m1
            Ne_DT = state_result[0]
    # return result_Pc_m1, Ne_m1, power_consumption_fan_m1, net_power_m1
    return [
                    round(power_consumption_fan_m1, 2),
                    round(net_power_m1, 2),
                    round(Ne_m1, 2),
                    round(m1, 2),  # 假设 data1[11] 是风机转速
                    round(Ta, 2),  # 假设 data1[10] 是环境温度
                    round(result_Pc_m1, 2),  # 转换单位为 kPa
                    pb_list,
                    ne_list,
                    state_result
    ]

if __name__ == "__main__":
    #ita_current, kc_current = 0.8057840192077641, 1.7175000000000007
    ita_current, kc_current = 0.720057840192077641, 2.4275000000000007
    results = optimize_for_speed(n_opt_initial, simtime=25, sample_time=5, h=0.2, tolerance=0.0005,
                                 ita_current=ita_current, kc_current=kc_current)
    print(f'实际数据仿真到稳态后负荷{results[2]}，背压{results[5]}')
    print(f'仿真到稳态后最终状态：{results[8]}')
    save_path = "pics/"
    os.makedirs(save_path, exist_ok=True)

    plt.figure(figsize=(12, 8))
    plt.title('实际背压仿真到稳态', fontsize=16)

    tick = [i for i in range(0, len(results[6]), 1)]
    # 绘制风机功耗增量
    plt.plot(tick, results[6], marker='o', linestyle='-', color='purple', )

    # 设置坐标轴标签
    plt.xlabel('仿真步数', fontsize=14)
    plt.ylabel('背压', fontsize=14)
    plt.xticks(fontsize=12)
    plt.yticks(fontsize=12)

    # 显示图形
    plt.tight_layout()
    plt.savefig(f"{save_path}/实际背压仿真到稳态.png", dpi=300)
    plt.close()

    # TODO：最优风机转速下，负荷仿真到稳态，排汽焓的变化，等等中间变量的变化
    plt.figure(figsize=(12, 8))
    plt.title('实际背压仿真到稳态过程中负荷变化', fontsize=16)

    tick = [i for i in range(0, len(results[7]), 1)]
    # 绘制风机功耗增量
    plt.plot(tick, results[7], marker='o', linestyle='-', color='purple', )

    # 设置坐标轴标签
    plt.xlabel('仿真步数', fontsize=14)
    plt.ylabel('负荷', fontsize=14)
    plt.xticks(fontsize=12)
    plt.yticks(fontsize=12)

    # 显示图形
    plt.tight_layout()
    plt.savefig(f"{save_path}/实际背压仿真到稳态过程中负荷变化.png", dpi=300)
    plt.close()
